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ABSTRACT 

Rapidly oscillating Ap stars exhibit an astrophysically interesting combination of strong, dipolar-like 
magnetic fields and high-overtone p-mode pulsations similar to the Sun. Recent time-resolved spec- 
troscopy of these stars unravelled a complex picture of propagating magneto- acoustic pulsation waves, 
with amplitude and phase strongly changing as a function of atmospheric height. To interpret these 
observations and gain a new insight into the atmospheric dynamics of roAp stars we have carried out 
2-D time-dependent, non-linear magneto- hydrodynamical simulations of waves for a realistic atmo- 
spheric stratification of a cool Ap star. We explore a grid of simulations in a wide parameter space, 
treating oscillations of the velocity, magnetic field and thermodynamic quantities in a self-consistent 
manner. Our simulations foster a new understanding of the infiuence of the atmosphere and the mag- 
netic field on the propagation and refiection properties of magneto-acoustic waves, formation of node 
surfaces, and relative variation of different quantities. Our simulations reproduce all main features of 
the observed pulsational behavior of roAp stars. We show, for the first time, that the overall depth 
dependence of the pulsations in roAp atmospheres is strongly infiuenced by the density inversion at 
the photospheric base. 

Subject headings: MHD - stars: magnetic fields - stars: oscillations - stars: chemically peculiar 



1. INTRODUCTION 
1.1. General properties of roAp stars 

Rapidly oscillating Ap (roAp) stars are main sequence, 
late-A chemically peculiar stars showing high-overtone 
p-m.o de pulsatio ns with periods between 6 and 22 min- 
utes (jKurtz Martinezll2QQQl : IKochukhovll2QQ8l) . Effec- 
tive temperatures of these stars span a relatively narrow 
range from 6400 to 8100 K. Their atmospheres exhibit 
large abundance anomalies, typical of the group of Sr- 
CrEu chemically peculiar stars. The light and iron-peak 
elements are normal or underabundant relative to the 
chemical composition of the solar photosphere, whereas 
concentration of the rare-e arth elements (REEs) is en- 
hanced by factors 10^-10^ (iRvabchiko va et al.ll200"i ). 

The presence of very strong, organized magnetic fields 
is another remarkable property of roAp stars. The 
mean field strength inferred from the polarization mea- 
surements and Zeeman splitting o f spectral lines usu- 
ally lies in the 1- 5 kC range (iMa thvs et al.' ' 19971 : 
iHubrig et all l2004bl : iKochukhov fc Ba gnulo 2006) but 
can reach 25-30 kG in exceptional cases (Kurtz et alj 
l2QQ6l : IFrevhammer et al] ^2008), making Ap stars the 
most extreme non-degenerate magnetic stars known. 
These fields have simple, dipolar-like topology, often 
significantly inclined with respect to t he stellar ro- 
tatio n axis (|Land street fc Mathvsl l2000l : [Bagnulo et alj 
I2OO2I ). Unlike rapidly evolving, dynamo-generated mag- 
netic fields of the Sun and late-type stars, the fields of 
Ap stars are most likely stable remnants of the fields 
acquired in the pre-main sequence evolutionary phase 
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(|Mossl 120041 : iBraithwaite fc Nordlundl I2006D . Ap stars 
also rota te significantly slower than their normal coun- 
terparts (|Abt fc Morrelllfl995l ). probably due to an en- 
hanced braking provided by a magnetized wind and a 
more efficient coupling with the circumstellar disk in the 
pre-main sequence stage (Stepieh 2000). 

The unusual atmospheric chemistry of Ap stars is 
attributed to the effects of atomic diffusion (|Michaudl 
llQTOT ). This slow process of chemical separation un- 
der the inffuence of competing forces of the radiative 
pressure and gravity is facilitated by the slow rota- 
tion and stabilizing effect of the strong magnetic field. 
Different diffusion velocities in the regions permeated 
by the horizontal and vertical magnetic field lines lead 
to substantial variation of chemical a bundances over 
the surfaces of magn etic A stars (Mich aud et aLlflOSlI : 
lAlecian fc Stiftll2004l ). These chemical starspots are re- 
sponsible for the spectacular rotational modulation of 
line profiles in the spectra of roAp stars, synchronized 
with the variation of the overall stellar brightness, photo- 
metric colors, magnetic field strength and the mean line 
of sight magnetic field comp onent (|Rvabchikova et all 
I1997I : IKochukhov et al.ll2004ar ). 

Chemical stratification in the line-forming atmo- 
spheric regions is another important consequence of the 
atomic diffusio n in Ap stars. Detailed empirical line 
profil e studies ([Bagnulo et al.l 1200 ll : iRvabchikova et al] 
I2OO2I ) and self- consistent theoretical diffusion calcula- 
tions ([Leblanc et al.|[2QQ9l ) show that the light and iron- 
peak elements are concentrated in the lower atmosphere 
layers, typically below continuum optical depth at A = 
5000 A logrs = -1 : 0. On the other hand, REEs, 
such as Pr and Nd which show many strong lines in 
the roAp spectra, are pushed by the radiative pressure 
to the optically thin layers above logrs = — 3 : — 4 
(Mashonkina et al. 2005, 2009). 

According to the current consensus, roAp pulsations 
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are driven by the n mechanism operating in the H I 
partial ionization zone with the additional influence of 
chemical gradients and ma gnetic quench ing of convec- 
tion (Balmfor thet alJl2QQll : lTheado et aLl [2009). The ef- 
ficiency of this excitation mechanism sensitively depends 
on the relationship between the magnetic field, convec- 
tion, and atomic diffusion. The physics of these processes 
is not entirely understood. Consequently, despite giving 
robust and informative pulsation frequency predictions 
fe.g.. ICunha et all l2QQ3l : iGruberbauer et all 12008). the 
current roAp excitation theories are unable to account 
for the H-R diagram distribution of 40 known roAp stars 
and cannot explain their difference from apparently con- 
stant but otherwise very similar non-pulsating Ap stars. 

Shortly after discovery of the p-mode pulsations in Ap 
stars a close connection between oscillations and mag- 
netic field became obvious. Pulsational amplitude and 
phase were observed to follow a rotational modulation 
curve correlated with the field strength variation. The 
times of the magnetic and pulsation maxima generally 
coincide, while pulsation phase shows a characteristic tt 
radian jump when the st ellar magnetic equ ator is closest 
to the line of sight fe.g.. lKurtz et all 1 19941 ) This behav- 
ior has inspired the oblique pulsator model (|Kurtzl fT982). 
which attributes roAp pulsation to a non-radial, axisym- 
metric, low angular degree modes aligned with the axis 
of dip olar magnetic field. Theoretical calculations (|Saiol 
l2005f ) and indirect im aging of the horizo ntal pulsation ve- 
locity field structure (|Kochukhovll2004l ) have vindicated 
the oblique pulsator model but also revealed that the 
modes are substantially distorted by the magnetic field 
and cannot be described by a single spherical harmonic 
function. 

1.2. Spectroscopic pulsation signatures of roAp stars 

The advent of efficient, high-resolution spectrographs 
at large optical telescopes allowed extending observa- 
tional analysis of roAp pulsations from traditional high- 
speed, broad-band pho tometry (see a thorough review by 
iKurtz fc Martinez 2000) to radial velocity (RV) measure- 
ments and studies of t he individual line profile variabil- 
ity (see reviews by lKochukhovll2007l . l2008l : iKurtd 12008). 
These investigations revealed a number of unusual atmo- 
spheric pulsation phenomena, not present in any other 
type of pulsating star and not anticipated theoretically. 

The most prominent and unique characteristic of 
the rapid spectroscopic variability of roAp stars is 
a large diversity of the pulsation amplitu des and 
phases of different spectral l ines (Kana an fc Hat zes 1998; 
Kochukhov fc Rvabchikova'2001; Mkrtichian et al. 2003; 
Elkin et al. 2005; Kochukhov 2006; Ryabchikova et al. 
2007a). The narrow cores of the hydrogen Balmer lines 
and absorption features of the rare-earth ions usually ex- 
hibit RV amplitudes from a few hundred m s~^ to several 
kms~^. On the other hand, lines of Ca, Si and iron- 
peak elements vary with amplitudes below ~ 50 ms~^. 
This amplitude discrepancy is frequently accompanied 
by a substantial difference in the pulsation phases. In 
extreme cases the phase lags be tween RV curves of dif- 
ferent lines reach 7r-27r radians dM krtichian et al]l2003l : 
iRvabchikova et al.ll2007"al : ISachkov et al. 2008). 

Equally puzzling is the pulsational bisector variation 
in many roAp stars. Time-series measurements per- 
formed at different intensity levels in the profiles of 



strong REE lin es and Ha core yield systematical l y dif- 
ferent results (ISachkov et all l2004l : iKurtz et all l2005l : 
iRyabchikova et al.ll2007al )bl). tvpicallv showing a large in- 
crease of the pulsation amplitude towards the line wings 
and occasional change of the pulsation phase across line 
profiles. 

These spectroscopic pulsation signatures of roAp stars 
are difficult to reconcile with the oscillation picture ex- 
pected for a normal, chemical l y- hom ogeneous stellar at- 
mosphere. 'Rv abchikova et ahl (j_2002) suggested that pe- 
culiar pulsation properties of roAp stars are related to 
the presence of strong atmospheric chemical composi- 
tion gradients, for which there is unequivocal evidence 
from observations and theory alike. Chemical stratifica- 
tion is particularly extreme for the rare-earth ions, which 
levitate so high in the atmosphere that their formation 
heights exceed those of the hydrogen line cores. Pulsa- 
tion waves propagate outwards in this chemically seg- 
regated atmosphere, gradually increasing in amplitude 
and consecutively perturbing the layers probed by the 
Fe-peak element lines. Ha and REE lines. Thus, time- 
resolved spectroscopic observations of roAp stars provide 
a vertical cross-section of p-modes and allow detailed 
study of the propagation and transformation of magneto- 
acoustic waves over a large range of geometrical heights. 

Several time-resolved spectroscopic studies suggested 
that the phase-amplitude behavior of roAp pulsations 
show characteristics of standing waves (nearly constant 
phase) in deeper layers and running wav es (changing 
phase) higher in the upper atmosphere (e.g.. 'Kurtz e t al.l 
2003). A few roAp stars exhibit node-like surfaces 
where the pulsation amplitude drops below the detection 
thresh old and the pulsation phase undergoes a tt radian 
jump (jRvabchikova et al . I [2007al : ISachkov et al.l [2008). 

High-quality spectroscopic time-series obtained for the 
strongest REE lines (Kochukhov fc Ryabchikova 20Ql|) 
demonstrated that roAp pulsational line profile varia- 
tion is also highly unusual. Weak lines of the singly 
ionized Nd and Pr, which form in the lower part of the 
REE-enriched cloud, exhibit modulation expected for a 
low angular degree, non-radial pulsation. However, the 
strongest REE lines originating in the uppermost layers 
show asymmetric blue- to-red running waves, which are 
incompatible with the line profile variability pattern pro- 
duced by an y non-radial, horizontal p ulsational velocity 
fiuctuation. iKochukhov et all (|2007l ) successfully mod- 
eled this phenomeno n by a combination of t he RV and 
line width variation. IShibahashi et al.l (|2008 ) suggested 
an alternative explanation involving very high amplitude 
shock waves, unresolved by observations. 

The question of the magnetic field variation with pul- 
sation ptmse_Jia^been_a topic of intense recent de- 
bate. iHubrig et al.l (|2004al ) predicted magnetic variabil- 
ity of the order of 6B/B ^ 0.1 but could not reli- 
ably detect pulsational field changes in any of the six 
roAp stars they surveyed with l ow-resolution spectropo- 
larimetry. Leone fc Kurtz* (2003) claimed a detection of 
6B/B = 0.2 pulsational field modulation in the REE 
lines of the roAp sta r 7 Equ. However, thes e results 
were not confirmed bv IKochukhov et al.l (|2004cl ) who es- 
tablished an upper limit of 6B/B = 0.02 for the same 
star based on a more extensive observational material. 
Subsequent se arches of the pulsational v a riation of mag- 
netic field by IKochukhov et al.l (|2004bl ): ISavanov et al.l 



Magneto-acoustic waves in roAp stars 



3 



(|2QQ6[ ): iKochukhov fc Wadd ([2003) also yielded null re- 
sults. 

1.3. Overview of MHD wave calculations for the Sun 
and roAp stars 

To understand the complex physics of waves in 
magneto- atmospheres, numerical modeling of wave 
interaction with magnetic field has become a preferred 
approach in recent years, with the development of 
large computational facilities. It has an advantage over 
the analytical models, allowing to include magnetic 
field configurations of nearly arbitrary complexity 
and to study compound effects of several different 
physical agents, such as radiative losses or non-linear 
effects. For the Sun, where localized magnetic struc- 
tures (sunspots, fiux tubes) can be nearly spatially 
resolved in observations and studied in detail, the 
numerical modeling of waves in the photosphere and 
chromosphere of these structures has been per f ormed 
for several decades, starting from e.g. Shibat^ (|l983l ) 
and exponen tially increasing i n tim e durin g the last 
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without touching the literature on waves in the corona 
where the physical conditions are significantly differ- 
ent from those in the lower atmosphere). Important 
physical conclusions have been reached by these models. 
It has been demonstrated that the behavior of waves 
depends crucially on whether they propagate in the 
magnetically dominated region or in the gas pressure 
dominating region. The separation between these two 
regions is defined by the value of the plasma parameter 
P = Pg/Pm^ the ratio between the gas pressure and 
the magnetic pressure. In the particular case of waves, 
a more adequate parameter is the ratio between the 
squared characteristic wave speeds, the sound speed 
and the Alfven speed, c\/v\. These two parameters 
are related as {c\/v\)/l3 = 7/2 = 5/6 (for an ideal 
monoatomic gas). In the magnetized stellar atmosphere 
the ratio c^/v\ decreases exponentially with height 
as the density decreases. In the solar case a typical 
situation is to have cs of the same order of va at the 
spectral line formation heights in the photosphere or 
low chromosphere, depending on the field strength of 
the magnetic structure. 

In a two-dimensional situation, when only fast and 
slow MHD modes exist, the slow mode is mainly mag- 
netic in the region where cs > va and is mainly acous- 
tic in the case cs < va- The fast mode nature is 
the opposite. A detailed description of the fast and 
slow modes behavior, refraction and mode convers ion in 
the so l ar atmosphere can be found in e.g. Bogdan et al.l 
(|2003r ): IKhomenko ColladosI (|2006r ). The slow mode 
is mostly field-aligned in the sense that its group ve- 
locity is parallel to the local magnetic field direction. 
The fast mode can propagate across the field. To es- 
tablish the wave type in observations it is important 
to verify the height of the cs = va level with respect 



to the line formation level. Recently a lot of attention 
has bee n drawn to the phenomenon of mode transfor- 
mation JZhugzhda fc Dzha lilovlll982l 119841 : ICallvl 120061 : 
iCallv fc ^oossens 2008) taking place around the cs = va 
level. Applied to the Sun, mode transformation is be- 
lieved to be responsi ble for the wave powe r reduction 
observed in sunspots (|Callv fc Bogdanlll997l ) and, possi- 
bly, for the wave power redistributi on over active regions 
(jHanasogd [20081 : iKhomenkol 120081 ). Depending on the 
inclination of the field with respect to the wave propa- 
gation direction, the transformation can be complete or 
not. But in any case, new modes appear after the trans- 
formation and the wave behavior around the cs = va 
level can be extremely complex since different modes, 
propagating up and down, can interfere. In a three- 
dimensional situation, the wave behavior is even more 
mixed, and our overall understanding is far fro m being 
complete (however, see lCallv fc Goossensll200"8l ). 

Regarding the magnetized atmospheres of the roAp 
stars, it has been generally accepted that the mag- 
netic field makes important effects on osc i llatio n prop- 
erties (see the recent reviews by ICunhal 120071 : iKurtd 
|2008[ ) . A number of theoretical studies have been carried 
out to advance understanding of the infiuence of mag- 
netic field on oscillation eigenfrequencies, eigenfunctions, 
and on the structure of the multiplets and frequency 
spacing:s fc )ziembowski fc Goode 1996; Bigot et al. 2000; 
I Cunha fc GQUgh.2000: ,.Saio fc G autschv 2004: Saio 2005. : 
ICunhall2006l ). These theoretical models were able to pro- 
vide a qualitative explanation of the behavior of the fre- 
quency spa cings due to t he action of the magnetic field 
(|Cunha fc Gou gh 2000: iSaio fc Gau tschv 20 0j). Contin- 
uing the work of Cunha fc Goughl (|2000 ). CimS] ^20M) 
used a variational approach to calculate the magnetic 
perturbations to eigenfrequencies. She assumed that 
magnetic perturbations are introduced only in a thin 
layer around cs = va (which was called "magnetic 
boundary layer") changing the wave vector and intro- 
ducing additional phase shifts. Note that a similar 
picture was used to explain the behavior of the helio- 
seismic waves in solar spots, c allin g the mag netic ef- 
fects "sur face effects" (iBraunI 119971: iLindsev fc BraunI 
l2005iJ!bl; iMoradi fcCallvl l2008l: iMoradi et all l2009l : 
Khomenko^IelaU 120091). ISaio fc Gautschvl (|2004f ) and 
Saio (2005) addressed the same problem with a differ- 
ent approach, by using expansio ns of the eigenfunctions 
into s pherical harmonics (see also lDziembowski fc Good^ 
Il996l ). In both approaches, it was found that, depending 
on its strength, the magnetic field causes cyclic effects 
on the pulsation frequencies. Due to the mode conver- 
sion at the magnetic boundary layer a part of the wave 
energy escapes to the high-frequency running acoustic 
waves (slow low-/? magneto- acoustic modes) or to the 
slow high-/? magneto- acoustic mode waves in the inte- 
rior, producing mode energy losses that can become im- 
portant at some particular frequencies. 

As was mentioned above, the amplitudes and phases 
of the observed waves depend on the atmospheric height 
in a complex way. This behavior was larg e ly no t un- 
derstood until the rec ent studies by ICunhal ('2006') and 
ISousa fc Cunhal (|2008l ). Unlike the Sun, oscillations ob- 
served in roAp stars come from a magnetically domi- 
nated region in the most part of the stellar atmosphere. 
Only in the deepest observed layers the gas pressure and 
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Lorentz forces can be of the same order of magnitude 
for the lowes t field strengths encountered in roAp stars 
(|Cunhall2QQ7l ). Due to the direct influence of the Lorentz 
force, running acoustic waves in the atmosphere are sup- 
posed to become field-aligned and their acoustic cut-off 
frequency changes as a function of the magnetic field in- 
clination as uoc cos 6>, where ujr^ is a cut-off frequency in the 
absence of magnetic fi el d (iDziemb owski fc Go odd 119961 : 
iCunha fc Goughl [2QQQh . iSousa fc Cunha (200^onsid- 
ered an analytical quasi one-dimensional model of the 
radial high-frequency pulsations (above the cut-off fre- 
quency) in a polytropic interior model of roAp stars 
matched to an isothermal atmosphere. They found that 
in the interior running magnetic waves with progres- 
sively decreasing wavelength exist (slow high-/3 magneto- 
acoustic modes) together with nearly standing acous- 
tic waves with much larger wavelength (fast high-/? 
magneto- acoustic modes). In the atmosphere the sit- 
uation is the opposite: running acoustic waves with 
the displacement parallel to the magnetic field coexist 
with the nearly standing magnetic waves with displace- 
ment perpendicular to the magnetic field (slow and fast 
low- j3 magneto-acoustic modes, respectively). A super- 
position of the line of sight projections of both wave 
components non-trivially depends on the magnetic field 
parameters and can possibly produce a complex spec- 
troscopic pulsation behavior, in qualitat ive agreement 
with observations (|Sousa fc Cunhal[2QQ8h . For oscilla- 
tions with frequencies below the cut-off frequency, the 
losses to the running acoustic waves in the atmosphere 
are restricted to latitudes where the inclination of the 
dipolar magnetic field is sufficient to lower down the cut- 
off according to uOc cosO. For high-frequency oscillations, 
only a part of the energy is removed due to the run- 
ning acoustic waves in the atmosphere and they can be- 
come trapped, unlike the non-magnetic case, where high- 
frequency waves are always propagating. The losses to 
the high-frequency running acoustic waves in the atmo- 
sphere were found to be the largest near the poles of 
the magnetic dipole where the field is nearly vertical. 
The losses to the running magnetic waves in the inte- 
rior are the largest for fields inclined around 30 degrees. 
These results are similar to those from local helioseismol- 
ogy inode conve rsion studies in inclined magnetic fields 
by ICiOWl ([20061) and iSchunker fc Cidhl (120061 Thus, 
except for regions near the poles, the high-frequency 
w ave energy is nearly t rapped. Extrapolating the results 
of ISousa fc Cunhal (|2QQ8) , waves with frequencies below 
the cut-off frequency are expected to be trapped at and 
near the poles but able to convert into running acoustic 
slow-mode waves with increasing inclination, so that the 
mode energy losses would be maximum in an annulus 
close to the poles (but not exactly at the poles). In all 
the cases, trapped waves concentrate near the magnetic 
equator. Note, that earlier works tried to explain the 
high-frequency wave reflec tions by an ad hoc chromo- 
spheric temperature rise (Gautsc hv et al.M l998). which 
contradicted observations since no evidence of chromo- 
spheres was found for the cool Ap stars. 

The observational evidence for the lack of magnetic 
field variation ov er pu l sation cycle seems to be confirmed 
theoretically by Saio (2005) who found that the mag- 
netic field pulsations do not exceed the hmit of 6B/B = 
10~^. However, these theoretical calculations relied on 



an approximate analytical description of the atmospheric 
structure and did not extend to the low optical depths 
where strongly pulsating lines are formed. 

1.4. Goals of our study 

Despite this recent progress, more theoretical work is 
needed to reach the complete understanding of the roAp 
star's pulsation properties and, in particular, to explain 
the unusual spectroscopic variability of these stars. In 
this paper we use two-dimensional numerical simulations 
to study the properties of magneto-acoustic waves in the 
atmospheres of roAp stars. We take advantage of the nu - 
merical code developed bv lKhomenko fc ColladosI (|2006l ) 
to model the wave propagation in magneto- atmospheres. 
This code was previously applied to a variety of problems 
of solar physics, such as the wave energy transport to the 
chromosphere in small-scale magnetic flux tubes, or he- 
lioseismic wave pr opagation though the magnetic struc- 
ture of sunspots (|Khomenko et al.l l2008aL 120091) ^__Our 
strategy is similar to the work of lSousa fc Cunhal (|2008l ) . 



However, since the equations are solved fully numerically, 
we were able to relax several approximations made by 
these authors. We assume the waves to be excited by a 
radial pulsation mode below the visible surface in a strat- 
ified non-isothermal atmospheric model permeated by a 
constant inclined magnetic field. The code solves non- 
linear equations, allowing for development of shocks in 
our simulations and treating the propagation of magneto- 
acoustic waves in a realistic cool Ap-star atmosphere. 
This approach allows us to relate interesting pulsation 
properties with the formation heights of important diag- 
nostic lines, thereby bridging the gap between the com- 
plex observational picture of roAp pulsation and theo- 
retical models, previously limited to simplistic analytical 
atmospheric structures. 

Our goal is to create a more generalized picture and to 
study in detail the effects of the magnetic field strength 
and inclination, and the temperature stratification on 
the observed pulsation amplitudes and phases. We ana- 
lyze the behavior of individual magneto-acoustic modes 
at different latitude locations and at different frequen- 
cies, both below and above the atmospheric cut-off fre- 
quency. The signatures of these waves in the disc- 
integrated signal are also studied. Our analysis is not 
limited just to the velocity oscillations, but we also 
study self-consistently the temporal behavior of thermo- 
dynamic parameters and the of magnetic field vector. 

Our simulations are suitable for detailed spectral syn- 
thesis and a direct comparison to observations. How- 
ever, we defer this work to subsequent papers. Below we 
concentrate on the analysis of the physical properties of 
magneto-acoustic waves under different conditions and 
try to address several key questions: why such a diverse 
pulsation wave behavior is observed in the atmospheres 
of roAp stars, depending on the spectral line (i.e. the 
height in the atmosphere); which physical mechanisms 
determine the presence of the node surfaces, rapid growth 
of the amplitude and change of the pulsation phase with 
height, mixture of an upward and a downward wave prop- 
agation; what kind of temperature, density and magnetic 
field variations accompanies RV oscillations, etc. 

2. LOCAL MHD SIMULATIONS 
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Fig. 1. — Equilibrium model atmosphere, (a) Temperature, (b) Density, (c) Solid line: acoustic speed; dotted lines: Alfven speed for 
the values of the magnetic field marked with numbers, (d) Solid lines: cut-off period Tc = 27r/(c(;c cos -0), where ip is the magnetic field 
inclination and lJc = jg/2cs- Dotted lines: periods of the driver. In all panels zero height corresponds to logrs = 0. The optical depth 
scale is given by the upper axis. 



The global magnetic field of roAp stars can be 
approxim ated t o a fi rst order as a magneti c dipole 
^Bagnulo ilaD Il995l : iKochukhov et all l2QQ4af ). The 
characteristic wavelength of oscillations can be estimated 
as A ~ c^-T, where T is the pulsation period. For a typi- 
cal model atmosphere of roAp star cs is approximately a 
few tens of km s~^, depending on the atmospheric height, 
and T is of the order of ten minutes, making the wave- 
length to be of the order of 10^ Mm. This number is 
negligible compared to the stellar radius and to the char- 
acteristic scale of the variations of the magnetic field. 
Based on these considerations we take a local approach 
in the simulations, assuming that the magnetic field is 
locally homogeneous and inclined. 

Given that the observed pulsations are of low degree 
^, we assume that at the base of our simulation domain 
(below the visible surface) the oscillations are excited by 
a plane-parallel perturbation with mostly vertical dis- 
placement. For the dipolar field and modes with zero 
azimuthal order (m=0) the pulsations are axisymmetric. 
The simulated waves propagate in a plane defined by the 
magnetic field vector and the gravity vector. For the 
dipolar field the pulsations are axisymmetric and their 
properties depend on the magnetic latitude, but not on 
the azimuthal angle. Apart from this, the properties 
of oscillations in the horizontal direction are not con- 
st rained by any assumpti on. 

Dziembowski fc Goodd (119961) : iCunha fc Gmighl 



([2QQQh : iCunhal ([20061 ): iSousa fc Cunhal ([2QQ8h also 



solved a local problem in a plane-parallel atmosphere 
with a locally uniform magnetic field, allowing os- 
cillations to have the displacement vector in the 
two-dimensional plane. The difference between our 
approach and the approach of these authors is that our 
analysis is not limited to low-amplitude oscillations and 
allows arbitrary variations of the velocity in horizontal 
direction in contrast to quasi-radial pulsations consid- 
ered in the aforementioned papers. We focus entirely on 
the propagation of waves in the stellar atmosphere and 
a small portion of sub-surface layers and do not consider 
excitation of the pulsations and mode transformation in 
the stellar interior. 

2.1. Numerical method 

We use the numerical MHD simulation code that was 
developed during the last years at the Instituto de As- 
trofisica de Canarias (lAC) to study the wave propaga- 
tion in different mags:netic structures. The descrip t ion o f 
the code can be f ound in iKhomenko fc Collados (20061) ; 
iKhomenko et al.l (|2008bl ). The code solves the basic non- 
linear equations of the ideal MHD, written in conserva- 
tive form: 

^^'V-{pV)=0, (1) 



dt 



P9: 



(2) 
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(4) 



where I is the identity tensor and E is the total energy: 



E=\pV^ 



7-1 



B2 

Stt 



(5) 



Ah other symbols have their usual meaning. In these 
equations the gravity g is taken constant with a value cor- 
responding to that at the stellar surface. In the present 
version of the code, the radiative energy losses, Q, can 
be included according to the Newton cooling law. How- 
ever, we postpone the study of non-adiabatic effects to 
the future. 

The code solves the non-linear MHD equations for per- 
turbations, that are obtained by subtracting the equa- 
tions of initial magnetostatic equilibrium from Eqs. [TJH 
Note that this procedure is not equivalent to a usual 
linearization of the equations since we maintain all the 
non-linear terms in the equations for perturbations. In 
the simulations described above we use the equation of 
state of an ideal gas. 

The boundary condition issue is very important for 
wave simulations since artificial refiections at the bound- 
aries can distort the physical results and limit the dura- 
tion of the simulations. Our code makes use of a Perfectly 
Matching Layer (PML) boundary condition (Berenger 
IT994) that efficiently absorbs the incoming waves at the 
boundaries and prevent their return into the simulation 
domain, using relatively few grid points. In the simula- 
tions described below, the horizontal domain boundaries 
are periodic and the top boundary has a PML layer. At 
the bottom boundary of the simulation domain we drive 
sine waves with a given period by means of the exact an- 
alytical solution for the radial acous tic-gravity waves in 
an isothermal stratified atmosphere (|Mihalas fc MihalasI 
[^): 

Vz = Vo exp (^) sm{ujt - k^z) (6) 



p CO 
P CJ 



l^lexp (^^) sin(cjt 
^|exp (^^) sin(cjt- 



kzz + (I)r) (7) 



kzz + ^p) (8) 



where H = P/{dP/dz) is the pressure scale height and 
the amplitudes and relative phase shifts of the pressure 
and density perturbations are given by: 



\R\ 



1 



1 




(9) 

(10) 

(11) 
(12) 



Given the wave frequency cj, the vertical wave vector kz 
is found from the usual dispersion relation for acoustic- 
gravity waves kz = \/ [oS^ — uS^jcs- The wave vector kz 
is zero if the driving frequency is below the local cut-off 
frequency. 

As the atmosphere at the bottom boundary is not 
isothermal, the chosen form of driver is not an exact solu- 
tion for this layer. The vertical dependence in equations 
dlHH]) would be slightly different including the tempera- 
ture gradient. We have checked the implications of this 
form of driving on the solution developed in the physical 
domain. We found that varying artificially the vertical 
dependence in the equations dSHH]) till the extreme case 
of oscillations in the non-stratified atmosphere has neg- 
ligible infiuence on the solution in the physical domain. 
This gives us confidence that the numerical solution is 
not biased by this particular form of driving. 

2.2. Model atmosphere and simulation setup 

As unperturbed equilibrium model atmosphere we 
adopt a model computed with the help of the 
opacity-sampling; ni odel atmosphere code LLmodels 
CShul vak et aIll2QQ4f ). This plane-parallel, LTE model 
has effective temperature Tgff = 7750 K and gravitational 
acceleration log ^ = 4.0 (cgs units). The model takes 
into account abundance anomalies typical of cool mag- 
netic Ap stars and is part of the grid presented by 
iKochukhov fc ShulvakI (|2QQ8[ ). The dipolar magnetic 
field is force-free and, therefore, does not infiuence the 
hydrostatic equilibrium of the model. 

Some parameters of the model are displayed in Fig. [TJ 
It gives the temperature, density, characteristic speeds 
and cut-of period stratifications as a function of geomet- 
rical height and continuum optical depth logrs. Several 
features of the model should be emphasized for better un- 
derstanding the simulations. The temperature has a very 
strong gradient at the surface decreasing almost a factor 
of three in 0.5 Mm. This has consequences on the critical 
period Tc and on the acoustic speed 05- , which both fall 
abruptly at the atmospheric base. The density stratifi- 
cation has a bump around zero kilometers (logrs = 0), 
which appears due to the hydrogen ionization in these 
layers. Due to this bump and due to the temperature 
gradient, we can expect strong refiection of waves com- 
ing from the interior to the surface. As a consequence 
of the density stratification, the Alfven speed (Fig. [it) 
has a complex behavior showing a local minimum around 
zero kilometers. 

In the present paper we perform a grid of simulations 
made for different combinations of three parameters: lo- 
cal magnetic field strength, its inclination and driving 
period. The magnetic field strength grid has 5 values: 
^0 = 0.5, 1, 2, 3 and 7 kG. The inclination grid has 4 
values: ip = 0, 30, 45 and 60 degrees to the local verti- 
cal. And the period grid has 3 values: 360, 600 and 800 s. 
This makes in total 60 models. In this way we span the 
typical range of magnetic field strengths and pulsation 
periods observed in roAp stars. The initial amplitude 
of oscillations at the base of the simulation domain was 
Vo = 100 m s~^ for all simulations. Through this pa- 
per we take positive sign for velocity corresponding to a 
redshift (downflow). 

Note that we do not introduce any latitudinal depen- 
dence of the driving velocity. We use the same form 
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Fig. 2. — Height-time variations of the longitudinal (parallel to the field, top) and transversal (perpendicular to the field, bottom) 
velocities for Bq = 0.5 kG, Tq = 360 s and inclinations equal to tp = 0, 30, 45 and 60°. Yellow and white colors mean positive velocities 
and dark red and blue colors mean negative velocities. The velocity contours of |y| = (1, 0.5, 0.3) km s"-*^ are overplotted as solid lines. 
Zero height corresponds to the photospheric base (logrs = 0). Dotted lines marked with numbers are contours of constant c^/i;^. Height 
dependences of the sound speed cs (solid line) and Alfven speed va (dashed line) are plotted over the left bottom panel. 




of the driver (Eq. [6H8|) for all inclinations of the mag- 
netic field, i.e. at all latitudes of the star. In other 
words, we study the response of the stellar atmosphere 
to an acoustic wave perturbation introduced from be- 
low. By doing so we neglect the possible presence of 
other wave types at heights of the lower boundary of 
our simulation domain. For the range of the considered 
magnetic field strength, the low boundary is the region 
where the asymptotic solutions for different decoupled 
modes still can not be safely applied. The weight of the 
different modes at these heights depends on their prop- 
agation, refiection and transformation in the deeper lay- 
ers and is apriory unknown. Several theoretical studies 
(e.g. ISousa & Cunha 2008) show that the presence of 
magnetic boundary layer has strong implications on the 



latitudinal dependence of the amplitudes and phases of 
waves at the bottom of the atmosphere. However, we can 
not use the theoretical eigenfunctions calculated in such 
studies directly in application to our simulations, as the 
atmospheric model considered in them is different from 
ours. Giving these limitations, the reader has to keep 
in mind that the results of our calculations, as for the 
relative amplitude and phase distributions of waves at 
different latitudes (Sect. [3TLQ|) . are strictly valid only as- 
suming acoustic waves at the bottom of the atmosphere 
with the latitude-independent properties. 

Fig. [It gives the Alven speed for the Bq grid values. As 
we see from the figure, already for a field strength as weak 
as 1 kG, the Alfven speed exceeds the acoustic speed in 
the whole height range, making the atmosphere magnet- 
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Fig. 4. Same as Fig.[2]and Fig.[3]but for = 7 kG and Tq = 360 s. 
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Fig. 5.— Same as Fig.[2]and Fig.[3]but for So = 7 kG and To = 800 s. 
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ically dominated. Only for Bq = 0.5 kG the plasma is 
pressure dominated below —1 Mm. Thus, we expect to 
see some mode transformation in the corresponding sim- 
ulations. 

In the magnetically dominated atmosphere, the cut- 
off frequency changes with the magnetic field inclination 
angle as the acoustic wave propagation becomes field- 
aligned. Fig. [TJi illustrates acoustic cut-off periods for 
the inclination grid of our simulations. This figure also 
shows the driving periods Tq. In our model atmosphere 
waves with Tq = 360 s are propagating at all heights. 
However, waves with Tq = 800 s are propagating only in 
the interior and become evanescent at the surface. Thus, 
some refiection of the 800 s-period waves is expected at 
the surface. 

Through the rest of the paper, we describe the results 
of our simulations in terms of fast and slow magneto- 
acoustic-gravity waves (note that we use a more con- 



cise notation, i.e. "magneto-acoustic"). These waves 
are the two possible solutions of the MHD equations ([Ti- 
ll]), one with a larger propagation speed and one with 
a lower propagation speed. In an academic case of a 
homogeneous atmosphere that can be described by sin- 
gle values of Alfven speed va and the sound speed cs, 
the waves maintain their fast or slow properties globally 
in the whole domain. As we consider the case of non- 
homogeneous vertically stratified atmosphere, the prop- 
agation speeds of the two solutions change with height. 
In the regions in the atmosphere where va ^ cs or 
Va ^ Cs the two solutions are expected to recover their 
asymptotic behaviour identified with classical slow and 
fast waves of the homogeneous atmosphere. In the in- 
termediate situation when va < cs^ one can still iden- 
tify the fast mode as mostly acoustic and a slow mode 
as mostly magnetic. The opposite is true in the region 
where va > cs (see, e.g. lCallvll2001l : iBogdan et al.ll2Q03l : 
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iKhomenko fc ColladosI [2QQ6I , where such a classification 
is applied). The modes can be considered as "coupled" 
in the region where va ^ cs- The later term means that 
their properties become close. Even in this situation we 
preserve the notation "fast"' and "slow" , as the propaga- 
tion speed of one solution is always larger than the other 
one, except when va = cs strictly. 

3. RESULTS 

3.1. Velocity variations 

According to Fig. [TJ the simulated part of the atmo- 
sphere is in the magnetically dominated regime down to 
field strength of 1 kG. In this regime the natural reference 
system for oscillations is the one defined by the magnetic 
field. Below in this section we discuss oscillations of ve- 
locity parallel to the field and perpendicular to the field 
{V\ong and Vtran, respectively). This way the two, fast 
and slow, magnet o-acoustic modes can be distinguished 
(|Khomenko fc Co llados 2006). 

As mentioned above, the clear distinction between the 
two modes can be done only in the situation where 
Va ^ Cs or Va <C c^. In the va ^ cs limit, the slow 
mode is mostly an acoustic-gravity wave guided along 
the magnetic field lines and is visible in the longitudi- 
nal velocity and in the gas pressure variations. The fast 
mode is magnetic and is visible in the transversal ve- 
locity as well as in the variations of the magnetic field 
vector. In the va <C cs limit, the slow mode (now mag- 
netic) is better visible in the transverse velocity varia- 
tions and in the magnetic field vector variations. The 
slow waves in this re gion are quasi- Alfvenic in character 
(|Bogdan et al.ll2003( ), thus, the principal motion is trans- 
verse to the magnetic field and the propagation direction 
is along the magnetic field. While the fast mode in the 
Va ^ Cs limit is mostly acoustic and has no relation 
to the magnetic field. It can have both longitudinal and 
transverse velocity components accompanied by pressure 
fiuctuations. In the intermediate situation (va < cs or 
Va < Cs) the modes can not be so clearly distinguished, 
but still preserve their mostly acoustic or mostly mag- 
netic behaviour, depending on the region. Note that we 
keep the discussion in terms of Mong and Vtran projections 
even in this intermediate situation since it is convenient 
to keep the same notation over the whole atmosphere. It 
allows the analysis of the different mode behavior to be 
done in a clear way. 

Figs. [2]-[5] show Mong and Vtran velocities as a function 
of time and vertical coordinate. We take for illustra- 
tion four representative cases of simulations. In the first 
case shown in Fig. [21 the magnetic field strength is weak 
(^0 = 0-5 kG) and the wave period (Tq = 360 s) is well 
below the acoustic cut-off period. In this simulation, half 
way though the atmosphere the high frequency waves 
meet the layer where cs = va and we expect to observe 
wave mode transformation and refiections. In the sec- 
ond case, shown in Fig. [3l the magnetic field is similarly 
weak, but the wave period is above the acoustic cut-off 
period at the upper part of the atmosphere (Tq = 800 
s), giving the possibility to study the behavior of the 
evanescent waves. The third case (Fig. H]) illustrates the 
propagation of the high-frequency waves (Tq = 360 s) in 
a strong magnetic field (^o = 7 kG). Finally, the last case 
(Fig. [5]) is for low-frequency evanescent waves (Tq = 800 
s) in the strong magnetic field. In all the cases we show 



simulations at four magnetic latitude locations where the 
inclination of the magnetic field to the local vertical is 
ip = (magnetic pole), 30, 45 and 60°. Note that the 
format of these figures is such that larger inclination of 
the ridges means lower propagation speeds of waves and 
that vertical ridges mean infinite propagation speed, i.e. 
stationary waves. 

High-frequency waves in the weak field 

Fig. [2^ shows that, as expected, in the vertical field 
{ip = 0°) only one type of wave exists with the velocity 
aligned to the field. The vertical driver excites vertically 
propagating fast magneto-acoustic mode that is essen- 
tially acoustic at the bottom of the simulation domain 
where cs > va- Propagating upwards this wave encoun- 
ters the layer where cs = va and is completely trans- 
formed into the slow magneto-acoustic wave (also acous- 
tic). This mode only has longitudinal velocity compo- 
nent. It can be appreciated how the propagation speed 
of waves decreases above km, in agreement with the 
rapid decrease of the acoustic speed (Fig. [2^). Fig. [6^ 
and e show the amplitudes and phases of the Mong and 
Vtran vclocitics for the simulations with Bq = 0.5 kG and 
To = 360 s. The amplitude of the longitudinal velocity in 
the ip = 0° case (black lines) does not increase monoton- 
ically, but has two local minima and a local maximum at 
height around —2 Mm. This structure is formed by par- 
tial refiection of the acoustic wave because of the strong 
density, temperature and, thus, sound speed, gradient at 
Z = km. As the frequency of the wave is above the 
cut-off frequency in the whole domain (see Fig. [Tj) this 
refiection can not be due to the cut-off frequency effects. 
An instructive discussion on the partial refiection of high- 
frequency (above the cut-off) wave s in the atmosphere 
of Sun and stars is presented in Balmforth fc Go ughl 
(1990). In this paper the authors consider, among others, 
the case of a polytropic interior matched to an isother- 
mal atmosphere with a discontinuous jump in the sound 
speed. They come to the conclusion that refiection of 
high-frequency waves in such a jump can be significant. 
Their schematic model gives a very good approximation 
to our model atmosphere where the jump in the sound 
speed is produced at the photospheric base due to the hy- 
drogen ionization. In the solar case, partial refiections of 
waves from the photospheric temperature gradient were 
also studied by e.g. Marmolino et al. (1993). After the 
refiection the waves become partially trapped between 
—3 and Mm. This is evident from their phase behavior 
(black line in Fig. [6^) that is almost constant at these 
heights. Below and above the trapped zone, the acoustic 
waves are propagating upwards and their phase is mono- 
tonically increasing with height. 

The case of the inclined field is substantially more com- 
plex. The vertical driver in the inclined field excites a 
mixture of the fast and slow magneto-acoustic modes and 
their relative contribution depends on the field inclina- 
tion. On their way up the fast and slow modes encounter 
first the transformation cs = va layer. In the case of 
the inclined field the fast-to-slow (and vice versa) mode 
transformation is not complete and depends on the wave 
frequency and the angle between the wave vector and 
the magnetic field vector. According to Callvl (|2006l ) the 
higher the frequency, the smaller is the cone of angles 
where the fast-to-slow mode transformation is effective. 
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The longitudinal velocity variations show that the fast 
mode generated by the driver in the cs > va region is 
partially transformed into the slow mode in the cs < va 
region (Fig. [2)3, c and d). The comparison of the veloc- 
ity amplitudes in Fig. [B] for the case of three different 
inclinations indicates that the part of the driver en- 
ergy that finally reaches the upper layer in the form of 
the slow mode might decrease with the inclination angle. 
This happens because fast-to-slow mode transformation 
at the Cs = va layer becomes less and less effective with 
increasing m ags:netic field inclination, in agreement with 
ICallyl (|2QQ6l ). Note, that, in principle, this general trend 
might be modified con sidering more per i ods in the sim- 
ulations. As shown bv lCunha fc Gougt] (|2QQQl ). the effi- 
ciency of the mode coupling and its relation to the incli- 
nation of the magnetic field depends on frequency in a 
highly non-linear manner. 

The transversal velocity variations in Fig. [2f, g and h 
show that the slow (magnetic) mode is produced by the 
driver in the cs > va region. Preceding the results from 
Sect 1321 the presence of the slow magnetic mode is not 
only deduced from the transversal velocity variations but 
also from the significant magnetic field variations present 
in this region (see Fig. [12] for Tq = 600 sec, the results 
are similar for Tq = 360 sec). The slow magnetic mode is 
partially transformed into the fast magnetic mode in the 
Cs < Va atmosphere. Due to the drastic increase of the 
Alfven velocity in the upper layers the propagation speed 
of the fast mode becomes very large. The fast mode is 
partially reflected from such strong Alfven speed gradi- 
ent region and passes again through the transformation 
region. This generates additional fast and slow modes 
that interfere with the upcoming waves. The downward 
slow mode propagation below —0.5 Mm is clearly seen in 
Fig.[2f and Fig. [6^ (red dashed line). The node surface is 
formed due to the interference of the upwardly propagat- 
ing slow waves and downwardly propagating slow waves 
slightly above the cs = va height. The exact position 
of this surface depends on the vertical wavelength of the 
waves. Depending on the field inclination, the mixture 
of the modes generated by the driver and propagating 
upwards, those produced after the multiple mode trans- 
formations and those partially refiected due to the strong 
temperature gradient creates complex interference pic- 
tures, particularly evident in the case of the largest field 
inclination tj; = 60° (Fig. [2]i). The number and loca- 
tion of the node layers depends on the field inclination, 
both for Mong and Vtran, and is particularly complex in 
the case of weak field and low period (and, thus, low 
wavelength) waves (Fig. [6]). 

Low-frequency waves in the weak field 

The case of the low-frequency waves propagation in 
the weak field is illustrated in Fig. [3] and Fig. [71 (panels 
a and e). In this case, apart from the transformation 
layer, we have a layer where the wave frequency becomes 
larger that the cut-off atmospheric frequency. For waves 
with To = 800 s, this layer is located in the atmosphere 
where cs < va and, thus, its location depends on the 
field inclination. 

In the vertical field (?/; = 0°, Fig. [Sfci) the propagation 
properties are similar to the case described above. The 
waves with 800 s period are refiected around Z = Mm 
due to the cut-off frequency effects. Due to this reflec- 



tion, a node structure appears at about —2.5 Mm and 
a standing wave pattern is formed by the interference 
of the upcoming and downward reflected waves. Above 
Z ^ Mm, the 800 s period waves are evanescent. The 
evanescent wave pattern is evident from the phase and 
amplitude behavior of waves in Fig.[7K and e. The phase 
is almost constant with height, except for the tt radian 
jump at Z ^ —2.5 Mm, and the amplitude increase is 
much less that in the case of 360 s period waves. One 
may note a propagating behavior of the waves above 2- 
2.5 Mm. 

In the inclined fleld, similar to the previous case, a 
mixture of the fast and slow magneto-acoustic modes is 
generated at the lower boundary of the simulation do- 
main, propagating upwards. After the transformation 
layer, these modes are partially transformed and par- 
tially transmitted. The slow magneto-acoustic modes in 
the Cs < Va atmosphere are now affected by the cut-off 
frequency effects. For the fleld inclination ofip = 30° the 
position of the cut-off height has been shifted to Z ^ 1 
Mm and the waves become evanescent above this height, 
unlike the case of Tq = 360 s. With increasing inclina- 
tion, the cut-off height is moved higher up in the atmo- 
sphere outside the simulation domain and the waves be- 
come propagating upwards. This behavior is clear from 
the phase variations above Z = 1 Mm (Fig. [7^). The 
change of the slope of the ridges in Fig. [3)3, c and d with 
increasing inclination is due to the projection effects of 
the fleld-aligned propagation of slow mode waves. 

Fig. [3f, g and h show how the slow magneto- acoustic 
waves^ generated by the driver at the low boundary in 
the Cs > Va region are partially transformed to the fast 
magneto-acoustic waves in the cs < va atmosphere and 
again adopt a quasi-evanescent behavior due to the rapid 
increase of the Alfven speed. After the partial reflection 
of the fast magneto- acoustic waves, downward propagat- 
ing waves are clearly seen below ca = va height. Sim- 
ilar to the To = 360 s case, these downward propagat- 
ing waves are slow modes generated after the secondary 
transformation at the ca = va layer. The phase and am- 
plitude behavior is qualitatively similar for the Tq = 360 
and 800 s waves (compare dashed curves in Fig.[6K, e and 
Fig.[7K, e). The different pattern and node locations are 
due to the different vertical wavelengths of these waves. 
The height pattern is much smooth and continuous for 
the waves with larger periods. 

High-frequency waves in the strong field 

The results of the simulations of the high-frequency 
wave propagation in the strong magnetic fleld are illus- 
trated in Fig. [Hand Fig. [6] (panels d and h). The case of 
the strong fleld is much easier to analyze since the wave 
propagation always occur in the magnetically dominated 
regime and there is no mode transformation. In the 
high-frequency case there are also no cut-off frequency 
effects and the slow (acoustic) waves are propagating at 
all heights. 

In the vertical fleld the driver generates slow magneto- 
acoustic waves. In the magnetically dominated atmo- 
sphere these waves are essentially acoustic and their 

^ Again, the slow mode in this region can be considered as mainly 
magnetic as the variations in the transversal velocity are accompa- 
nied by significant magnetic field variations. 
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Fig. 6. — Amplitudes (top) and phases (bottom) of the longitudinal (solid lines) and transversal (dashed lines) velocities as a function 
of height for simulations with Bq = 0.5, 1, 3 and 7 kG and at different inclinations ip (indicated in the figure). The results are for the 
oscillation period To = 360 s. The amplitude curves for each ip are separated by 0.5 km/sec for better visualization. The function fitted to 
the simulations is Asin(c<;t + V^)? so that larger phase ip indicates a later time of pulsation maximum. 



propagation properties are identical to the case of Bq = 
0.5 kG, as the vertical magnetic field does not affect to 
the first order the properties of the vertically propagat- 
ing acoustic waves (see, however, Roberts 2006). In the 
inclined field, in the absence of the mode transforma- 
tion, the amplitudes of the slow modes that reach the 
upper boundary are very similar. The propagation of 
these waves is only affected by the strong gradient in the 
sound speed, producing partial refiection around Z = 
Mm. As the wave propagation is field-aligned, they at- 
tack the layer of the strong sound speed gradient with 
a different angle. We are not aware of any theoretical 
developments of the wave refiection on the sound speed 
gradient in the presence of the magnetic field. Based on 
the analogy with rays in geometrical optics (that applies 
in the case of high-frequency short spatial wavelength 
waves) we can speculate that the angle of refiection may 
depend on the angle of incidence to the refiection layer. 
This would explain different interference pattern formed 
by the upcoming and refiected slow waves for the field 
with different inclinations (compare Fig. [4^ to d). 

Fig. [4f, g and h show the fast magneto-acoustic modes 
generated by the vertical driver. The relative contribu- 
tion of the fast mode gets larger with the field inclina- 
tion t/j. Unlike the case of the weak magnetic field, the 
phase and amplitude of this mode remain constant with 
height. This happens because the propagation speed of 
this mode is extremely large and its vertical wavelength 



is much larger than the size of the simulation domain. 

Low-frequency waves in the strong field 

Finally, Fig. [5] and Fig. [71 (panels d and h) illustrate 
the propagation of the low-frequency waves in the strong 
magnetic field. This case is qualitatively similar to the 
high-frequency wave propagation in the strong field, ex- 
cept for the presence of the cut-off layer and a much 
larger vertical wavelength of waves. 

The slow magneto-acoustic waves generated by the 
driver are now refiected due to the cut-off. This hap- 
pens at different heights, depending on the field inclina- 
tion. As follows from the phase dependence on height in 
Fig. [7h, a standing wave pattern is formed in the simula- 
tions with =0, 30 and 45°. Unlike this, the 60° inclined 
field already decreases the cut-off frequency sufficiently 
enough to allow for the wave propagation. The phase of 
the slow-mode waves in the i/j = 60° case continuously 
increases with height. Note also, that for the same rea- 
son, their relative amplitude increase is more significant 
that in the case of the slow-mode wave propagation in the 
less inclined magnetic field. The height of the node layer 
increases with increasing field inclination. This happens 
because the waves are refiected at progressively larger 
heights. 

The fast magneto-acoustic waves generated by the 
driver are identical to the case of Tq = 360 s. 




Comparative study of the amplitudes and phases 

Figs. [HI [71 and [8] show the height dependence of the 
velocity amphtudes and phases for the different field 
strengths and inclinations. The first two figures were 
partially discussed above. Please recall that, as explained 
above in Sect. 12. 2[ the comparison of the amplitudes and 
phases should be done keeping in mind the particular 
form of the driver at the lower boundary (i.e. acoustic 
wave with properties independent on the magnetic field 
inclination). The relation between the amplitudes of the 
fast and slow modes for a given inclination ip^ as well 
as the relative amplitudes of waves at different ip can be 
modified assuming another form of the driver. 

Fig. [6] shows that for the frequencies well above the 
cut-off frequency (Tq = 360 s) the field-aligned velocity 
variations have a character of the upward propagating 
waves for the field strengths above 1 kG. Only for the 
lowest field strengths in our simulations, Bq =0.5 kG, the 
phase behavior with height is more complex and down- 
ward wave propagation is observed for some angles. This 
complex behavior is produced by the wave transforma- 
tion and mode mixing at the cs = va layer. It means 
that for the weakest fields the high-frequency waves can 
be trapped, unlike the pure acoustic case. While the 
vertical propagation is unaltered by the presence of mag- 
netic field, the magnetic effects become more pronounced 
with increasing inclination. For fields above 1 kG the 
waves propagate in the atmosphere dominated by the 
magnetic field and their behavior is more simple. The 



velocity component perpendicular to the magnetic field 
shows the behavior characteristic of standing waves for 
the field strength above 3 kG. In general, both amplitude 
and phase curves get smoother and less node structure 
is formed with increasing field strength. 

Fig. [71 shows that the field-aligned velocity variations 
for frequencies close to the cut-off frequency (Tq = 800 s) 
mostly behave as evanescent trapped waves. The only 
exception is the situation of the strongly inclined field 
sufficient to lower down the cut-off frequency. According 
to Fig. [7J the low-frequency waves are propagating for 

= 60°. It means that the low-frequency waves that 
would be trapped in the absence of the magnetic field can 
escape at latitudes near the magnetic equator if the field 
strength is sufficiently high (above Bq = 0.5 kG in our 
simulations). Again, due to the mode transformation, 
the height dependence of the amplitudes and phases is 
more complex in the Bq = 0.5 kG case. The transversal 
velocity variations show downward wave propagation for 
the weakest fields and standing waves for the stronger 
fields. 

Velocity variations for the Tq =600 s period show an 
intermediate picture between the Tq =800 and 360 s 
cases. In the Bq = 0.5 kG field case multiple reflections 
and phase mixing are present. The cut-off period is just 
around 600 s at the very top of the atmosphere (Fig. [T]). 
Because of that the waves in the ip = 0° case are reflected 
and a standing wave pattern is formed. In all other cases 
the waves are mostly propagating upwards. The ampli- 
tude increase with height and the location of the nodes 
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is more similar to the 800 s case, as the wave frequency 
is closer to the cut-off frequency and the wavelength is 
much larger than in the To =360 s case. 

Fig. [8] give the comparison between the velocity am- 
plitudes at the top of the atmosphere as a function of 
the field strength for the different inclinations and peri- 
ods. It provides information of how much of the initial 
driver energy reaches the top layer and how much of it 
is retained below due to different processes. First we see 
that for the vertical field the amplitudes of the longitu- 
dinal velocity are field-independent. The absolute values 
of the amplitudes in the vertical field case decrease with 
increasing the period. Both properties are in agreement 
with the analytical models, as the vertically propagating 
acoustic waves are not affected by the presence of the 
vertical field. Is gives us an additional confidence in the 
validity of the simulations. 

For the inclined field, the Mong amplitudes depend on 
the field and this dependence is different for the differ- 
ent periods. For the high-frequency waves (Tq = 360 
s, left panel of Fig. [8|) the amplitudes decrease with the 
field inclination for all ^o- This effect is due to non- 
linear behaviour of waves. The velocity amplitude of 
the linear the high-frequency waves scales as exp(z/2i^) 
with height. Preceding the results from Fig.[TOl the non- 
linearities become important at lower heights for propa- 
gation at larger inclination angles. The saw-tooth profile 
develops and the amplitude increase with height becomes 
slower than exponential. Thus, in the high-frequency 
case, the largest amplitudes should be observed near the 
magnetic poles. This situation is the opposite to the low- 
frequency wave case (Tq = 800 s, right panel of Fig. [8]). 
The velocity amplitude for the linear waves with frequen- 
cies below the cut-off frequency scales with height as 
exp(z/2i^ — {ujI — uj'^)/cs)^ making its increase with 
height weaker than in the high-frequency case. However, 
as the effective cut-off frequency becomes lower for the 
waves propagating along the inclined field, the second 
term in the exponent becomes smaller and finally dis- 
appears for sufficiently large inclinations as the waves 
become propagating again. As follows from the right 
panel of Fig. [HI this effect is dominating over the effect 
of non-linearities. Thus, in the low- frequency case the 
longitudinal velocity amplitudes increase with ijj and are 
the largest near the equator of the magnetic dipole. The 
case of To = 600 s (middle panel of Fig. [8|) is the inter- 
mediate between these two. Note that for = 60°, the 
amplitudes are in fact the same for all Bq and Tq. 

For all periods, the amplitudes of the longitudinal ve- 
locity at the top of the atmosphere are the lowest for the 
weakest field strength Bq = 0.5 kG. At the same time, 
the amplitudes of the transversal velocity are the highest. 
This is a direct effect of the mode transformation, as the 
part of the energy of the fast mode cs > va waves excited 
at the lower boundary does not reach the top in the form 
of slow mode acoustic-like waves but it is transformed 
to fast mode magnetic waves observed in the transver- 
sal component of the velocity. Note that the effects of 
the mode transformation become more pronounced with 
increasing inclination, as the fast-to-slow mode transfor- 
mation becomes less effective. In the Bq = 0.5 kG case 
the dependence of the longitudinal velocity amplitudes 
on is qualitatively the same for different periods. With 



increasing field strength, the effects of the mode transfor- 
mation disappear and the amplitudes become distributed 
according to the acoustic cut-off frequency effects as de- 
scribed above. It means that, if the field is sufficiently 
weak and the mode transformation can happen in the 
atmosphere, the amplitudes of the longitudinal velocity 
are the largest at the pole of the magnetic dipole, both 
for the high-frequency and for the low- frequency waves. 
Note that in all cases the transversal velocity amplitudes 
at the top of the atmosphere are very small compared to 
the longitudinal velocity and would contribute very little 
to the disc-integrated signal. 

In summary, in our simulations, all the waves in the 
weak-field case are partially trapped near the magnetic 
equator. In the strong-field case, the low- frequency waves 
are trapped near the poles due to the acoustic cut-off 
effects, but are propagating near the magnetic equator. 
The high-frequency waves in the strong field case are not 
trapped at all magnetic latitudes. 

Finally, Figs. [9l and [TOl describe the non-linear behavior 
of waves in our simulations. Fig. [9] shows examples of 
the temporal variations of the longitudinal velocity at 
the top of the atmosphere logrs = —5.8 measured in 
units of the local sound speed for three different periods. 
The saw-tooth shape of the temporal profile is evident 
for the smallest period Tq = 360 s. The profile is still 
asymmetric for Tq = 600 s but is very close to the normal 
sine- wave shape for Tq = 800 s. The amplitudes vary 
from around 0.25 x cs for Tq = 800 s to 0.7 x cs for Tq = 
360 s, but stay all the time below the local speed of sound. 
According to a classical definition of a shock wave, the 
amplitude of the disturbance should be above the local 
sound speed. Thus, the saw-tooth waves observed in our 
simulations are not proper shocks, but simple non-linear 
disturbances produced after the wave has run a certain 
distance over the atmosphere. 

Fig. [To] gives the height of the saw-tooth wave forma- 
tion as a function of the magnetic field strength for dif- 
ferent inclinations and periods in the simulations. This 
height was defined in the way explained below. In Fig. [9] 
we marked three extreme points of the temporal pro- 
file A, B and C. If the slope between the points B 
and C is twice steeper than the slope between A and 
B we say there is a saw-tooth profile. The height of 
the saw-tooth profile formation depends on the wave pe- 
riod and its amplitude (Priest 1984; Mihalas & MihalasI 
!l984l ) and is typically larger for long-period waves. In 
our simulations with the initial amplitude of the veloc- 
ity at the base of the atmosphere of Vq = 100 ms~^ the 
saw-tooth profile formation does not occur below 2.2 Mm 
(logTs = —4.5). Fig. [To] shows that indeed, the height of 
the saw-tooth profile formation is the lowest for the lower 
period waves. The saw-tooth profile formation is less ef- 
fective for the weakest field as the wave transformation 
effects prevent the amplitudes from rising sufficiently to 
produce such profiles. In all the cases (except Bq = 0.5 
kG), the saw-tooth profile formation height is the low- 
est at high inclinations because the distance over the 
which the waves run to reach the same height is larger 
for the larger inclinations. 

3.2. Variations of thermodynamic parameters and 
magnetic field 
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Fig. 8. — Amplitudes measured at the top of the atmosphere at Z=2.8 Mm (corresponding to logrs = —5.2) as a function of the field 
strength for different inclinations, as indicated on the figure. Left: period To = 360 s; middle: period To = 600 s; right: period To = 800 s. 
Solid curves: Viong; dashed curves: Vtran- 
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Fig. 9. — Example of the temporal variations of the Viong at 
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^ = and three different periods To =360, 600 and 800 s. The 
amplitudes of Viong cire measured in units of the local sound speed 
cs- The three curves are separated by 1 x cs- 



Our simulations allow to calculate self-consistently 
the variations of the thermodynamic parameters and 
the magnetic field accompanying the velocity variations. 
Fig. inland Fig. [T2l illustrate the pressure, temperature, 
density and horizontal magnetic field variations for the 
simulations with 5o = 3 kG and Tq = 600 s. Unlike the 
case of the velocity data, we do not show simulations for 
the extreme cases of Bq and Tq but rather discuss the 
intermediate situation. The behavior of the thermody- 
namic variables and the magnetic field is rather similar 
in all the simulations with only some quantitative differ- 
ences. 

Fig. [TT] shows that the propagation of the pressure, 
temperature and density disturbances from the deep lay- 
ers to the surface is not continuous. There is an ampli- 
tude and phase jump around Z =0 Mm height, where 
the properties of the atmosphere change in an abrupt 
way. An interesting feature is observed in the density 
and temperature variations at this layer. There is a very 
localized amplitude increase (a "bump") present for all 
field inclinations. The nature of this localized enhance- 
ment is probably related to the density jump located 
around Z =0 Mm. Fig. [12] shows that the amplitudes 
measured relative to the static values are very high at the 
"bump" and can reach some 5-10 percent. The phase of 
perturbations at the "bump" seems to suffer an abrupt 



180° switch. The dependence of the amplitudes at the 
"bump" on the magnetic field strength 5o and inclina- 
tion il) is rather complex. For the weak field {Bq = 0.5 
kG) the amplitude of both 5T /T and 5p/ p increases with 
the field inclination. On the other hand, the opposite is 
true for larger fields {Bq > 3 kG). Note from Fig. [T2] 
that the "bump" is located in the region of formation of 
continuum intensity around logrs = 1 : — 1. Therefore, 
we can expect that these strong temperature variations 
will lead to continuum intensity variations which would 
dominate the total fiux changes of the star. The "bump" 
is absent in the relative pressure variations. 

The amplitude of variations of the thermodynamic 
parameters increases exponentially with height above 
logTs = — 1 and reaches some 15-20% at the top of the 
atmosphere. In all the cases the amplitudes are smaller 
for larger inclinations. As follows from the time- height 
pattern of the pressure and density variations for larger 
^l) (Fig. [TT]), at heights above the "bump" both the slow 
acoustic-like modes and fast magnetic modes contribute 
to these variations. The standing wave pattern char- 
acteristic for the fast modes is especially visible in the 
ip = 60° case at the first 20 min of the simulations. In- 
terestingly, the presence of the fast mode seems not to 
affect the temperature variations. Similar to the velocity 
oscillations, the oscillations of the thermodynamic pa- 
rameters are more affected by the magnetic field in the 
weak-field case (Fig. [121 l^ft column of panels). Several 
node structures are formed in the 5P/P and 5p/ p varia- 
tions for ip = 45° and 60°. With increasing field strength, 
the height variations of the amplitudes become very close 
for all field inclinations. 

The last row of panels in Fig. [TT] and Fig. [12] give the 
time-height plots and the amplitudes of the variations 
of the horizontal component of the magnetic field vector 
SBx . The variations of the magnetic field are not aware of 
the jump of thermodynamic parameters and their prop- 
agation is more continuous. There are no variations of 
the vertical magnetic field component. The amplitudes of 
variations of the horizontal magnetic field component are 
very small and do not exceed 3-4 G. The SBx amplitude 
decreases with increasing the magnetic field strength so 
that the magnetic field variations are almost absent in 
the Bq = 7 kG case. The 6Bx variations are larger in 
the deep layers and decrease abruptly above logrs = — 2. 
The amplitudes do not have a pronounced dependence 
on ip. The variations are absent for = 0° as the wave 
propagation is purely acoustic-like. 
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Fig. 10. — Height of the saw-tooth wave formation in Viong as a function of the field strength for different inchnations, as indicated 
on the figure) . The right vertical axis gives the optical depth logrs for the corresponding geometrical heights. Left: period To = 360 s; 
middle: period Tq = 600 s; right: period Tq = 800 s. 
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Fig. 11. — Height-time dependences of relative variations of the gas pressure (first row), temperature (second row), density (third row) 
and horizontal component of the magnetic field vector (last row) for the simulation with Bq = 3 kG and To = 600 s and inclinations equal 
'0 = 0, 30, 45 and 60°. Yellow and white colors mean positive values while dark red and blue colors mean negative values. The contours 
of \SP/P\, \ST/T\, \Sp/p\ = (4, 10, 20) % are plotted as solid hues over the first three panel rows. The contours of SBx = (0.2, 0.5, 1) G 
are plotted over the last row of panels. Zero height corresponds to the photospheric base at logrs =0. 



Finally, Fig. [13] illustrates variations of the gas den- 
sity, temperature and vertical velocity component at sev- 
eral different optical depths for the representative simu- 
lation set (To = 600 s, 5o = 3 kG, i/j = 0°). Unlike 



the previous height-time and amplitude- height plots, this 
figure shows absolute variations of the thermodynamic 
quantities. This format emphasizes strong, anti-phase 
changes of the temperature and density at the base of 
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Fig. 12. — Amplitudes of the relative variations of the gas pressure (first row); temperature (second row); density (third row) and 
horizontal component of the magnetic field vector (last row) as a function of the optical depth logrs. Columns from left to right are for the 
magnetic field strengths = 0.5, 1, 3 and 7 kG. Different curves are for different inclinations, the color coding is indicated in the figure. 



the photosphere. In the higher layers the amphtude of 
the temperature variations first rapidly decreases and 
then gradually increases again, reaching values compa- 
rable to the amplitude at logrs = only in the upper- 
most layers. The variation of the vertical velocity com- 
ponent rapidly grows above logrs — 1 and acquires a 
noticeable non-sinusoidal shape at logrs = — 5. In these 
highest layers the density variation appears to be rel- 
atively unimportant. In the upper atmosphere, above 
logr5 = — 1 the phase shifts between temperature, pres- 
sure and density variations follow the behavior character- 
istic f or acoustic-gravity waves dominated by buoyancy 
force (|Mihalas fc Miha'iaslll984i ). The temperature leads 
downward velocity by some 0.12-0.25 of pulsation pe- 



riod (45°-90°), depending on height, while the density 
lags behind the velocity by nearly the same amount. At 
layers around the "bump" at logrs = 0, the phase re- 
lations are switched by 180°. Of course, details of this 
pulsation picture depend somewhat on the assumed pul- 
sation period and magnetic field characteristics but the 
main features persist through all simulations. 

4. GLOBAL MODEL 

In addition to the extensive grid of local simulations for 
different values of the pulsation period, magnetic field 
strength and orientation, we performed a more restric- 
tive calculation aimed at reproducing the disc-integrated 
signal of the oblique non-radial pulsations observed from 
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Fig. 13. — Absolute variations of the vertical velocity, temper- 
ature and gas density for the simulations with Tq = 600 s and a 
vertical field with Bq = 3 kG. Pulsation curves for different opti- 
cal depths are shifted vertically. Positive velocity direction corre- 
sponds to redshift (downfiow). 



different aspect angles. As discussed in Sect. II. ![ the 
horizontal structure of the roAp pulsation generally re- 
sembles the axisymmetric £ = 1 pulsation mode, aligned 
with the dipolar magnetic field. Therefore, we adopted 
the field components: 

Br = cos (13) 
Be = 0.5Bd sin i9, 

where 6 is the latitude measured from the magnetic pole 
and 5d is the polar magnetic field strength. The strength 
and inclination of the magnetic field at latitude are 
given by: 

^0 = ^52+^ (14) 

i/j = arccos Br/Bo . 

We considered two dipolar field strength values, B^ = 
0.8 and 3 kG, and adopted the latitude-dependent driver 
amplitude, Vd = Vo cos6>, according to the simplest form 
of the oblique pulsator model. Apart from that we did 
not introduce any latitudinal dependence in the proper- 
ties of the driver as given by the equations (|6H8]). As 
explained in Sect. 12.21 the results of the disc integration 
are dependent on this particular form of the driver. This 
should be kept in mind analyzing the disk-integrated 
variations presented below. Similar to the local models, 
we used Vq = 100 m s~^. The driving period was set to 
600 s. Simulations were performed for an equidistant lat- 
itude grid in cos6>, varying from cos 6 = 0.95 to 0.05. For 
the purpose of disk integration, we divided each hemi- 
sphere of the star into 10 latitude zones, corresponding 
to this 6>-grid. Each latitude was further partitioned into 
60-100 elements, giving us a total of 1148 surface zones. 

Four angles are required to fully describe the ori- 
entation of the stellar rotation and magnetic axes 



with respect to the obser ver (e.g., see Fig. 6 in 
iPiskunov fc Kochukhovl[2002l ). However, here we are not 
interested in the transverse magnetic field component. 
Moreover, we limit ourselves to the investigation of the 
disk average of the line of sight velocity component (ra- 
dial velocity), which is not altered by the st ationary ve- 
locity field produced by the stellar rotation (jKochukhovl 
l2005f ). With these simplifications the only relevant pa- 
rameter is the angle 7 between the pulsation axis and 
the line of sight. 

We performed the disk integration by transforming the 
surface pulsation structure from the stellar to the ob- 
server's reference frame, rotated the coordinate system 
according to the angle 7 and summed, with an appro- 
priate weight, different pulsation quantities over all vis- 
ible surface elements. This procedure was applied in- 
dependently for each depth layer and simulation time 
step. The weight function took into account the pro- 
jected areas of the surface zones and a linear limb dark- 
ening with the coefficient u = 0.5. We have verified that 
the line of sight magnetic field component, Aos, obtained 
with this numerical scheme reproduces the longitudinal 
field calculated w ith the well-known analytical formula 
(jLerov et al.lll994f ) to within 3%. 

Fig. [m presents the height-time variations of differ- 
ent quantities for the global model with B^ = 3 kG. In 
this plot a positive radial velocity T^ios corresponds to 
the redshift as seen by the observer. On the other hand, 
a positive longitudinal field Bio^ points towards the ob- 
server, according to the standard definition used in solar 
and stellar magnetometry. 

As expected, the amplitude of pulsations depends on 
7. However, the overall depth dependence is dominated 
by the pulsations close to the magnetic pole and hence 
is similar for different angles between the line of sight 
and the pulsation axis. The radial velocity variations 
are qualitatively similar to the behavior of Mong for the 
low ?/^, strong-field, low- frequency local models (Fig. [5]). 
Thus, the radial velocity pulsations are dominated by the 
slow acoustic-like modes in all the visible atmospheric 
layers. This result is easy to understand giving the fact 
that the amplitudes of the transversal velocity variations 
produced by the fast magneto-acoustic modes in our sim- 
ulations are always smaller than the amplitudes of the 
longitudinal velocity variations due to the slow magneto- 
acoustic modes, except for the magnetic latitudes with a 
very inclined magnetic field i/j = 40°— 60° (see e.g. Fig.[7j). 
The node structure present in the radial velocity oscil- 
lations at logTs = is the remnant of the nodes due to 
the wave refiection on the density jump at the photo- 
spheric base. This node is persistent at all latitudes and 
does not disappear in the integrated velocity data. Apart 
from this node, the wave amplitude in the high layers in- 
creases rapidly with height reaching around 1 km s~^ at 
the top of our atmosphere. 

Fig. [m shows that the prominent temperature fiuc- 
tuations around logrs = are retained in the disk- 
integrated picture. Thus, this pulsational feature will 
have an important infiuence on the disk- averaged observ- 
ables. The amplitudes of the temperature fiuctuations 
reach a few percent, a value similar to individual vari- 
ations at each particular magnetic latitude. The tem- 
perature variations are the highest when the pulsation 
axis and the observational line of sight are co-aligned. 
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Fig. 14. — Height-time variations of the radial velocity (first row), relative temperature (second row), relative pressure (third row) and 
the line of sight magnetic field component (bottom row) for the global model with = 3 kG and To = 600 s. Columns show results for 
different values of the angle 7 between the pulsation axis and the line of sight. The contours of |Viosl = (0.05, 0.1, 0.2, 0.4, 0.8) km s"-*^, 
\ST/T\ = \SP/P\ = (1, 3, 5, 8) %, l^iosl = (0.05, 0.1, 0.2, 0.4) G are plotted as solid lines. The horizontal dotted lines correspond to the 
optical depths logrs = 1.5, 0, —2 and —4. 



The variations of the mean field modulus (not shown) 
exhibit a height dependence comparable to Bio^ and an 
amplitude below 1 G. 

The global model results obtained for the dipolar field 
with 5d = 0.8 kG are similar to those described above. 
The main difference is a higher amplitude of the relative 
temperature variations at the photospheric base (up to 
5%), an increase of the magnetic field fiuctuations by a 
factor of two and a steeper change of the pulsation veloc- 
ity phase with depth, reminiscent of the Mong behavior 
in Fig. El 

5. DISCUSSION AND CONCLUSIONS 

In this paper we have described the results of the sim- 
ulations of magneto-acoustic wave propagation in the at- 
mospheres of roAp stars for a wide grid of magnetic field 
strengths and wave periods. Our simulations allowed us 
to obtain a more generalized picture of the wave behav- 
ior, refiection, mode transformation, formation of evanes- 
cent waves and node layers. Our main findings can be 
summarized in a following way: 

• The atmospheric regions above logrs = of the 
considered roAp star model are magnetically dom- 
inated already for the field strengths as low as 0.5 
kG. Both fast and slow magneto- acoustic modes 
are present simultaneously in these layers. Ex- 
cited by vertical driving in our simulations, the slow 



mode waves dominate in the regions of nearly ver- 
tical field close to the magnetic poles. The fast 
mode waves only have comparable amplitudes to 
the slow mode waves at latitudes close to the mag- 
netic equator, corresponding to larger that 50°- 
60° (depending on the frequency). 

• The slow mode waves are mostly propagating up- 
wards and are field-aligned and the fast mode waves 
acquire a standing wave pattern in the upper at- 
mosphere due to their very large vertical wave- 
length produced by the rapid increase of the Alfven 
speed. This smooth behavior is disturbed by differ- 
ent refiection layers: cut-off layer, cs = va trans- 
formation layer and the density "bump" around 
logTs = 0. Due to refiections on (or above) these 
layers, trapped waves or waves propagating down 
are produced in the regions below the visible sur- 
face in our simulations. 

• The node layers can form due to the wave refiec- 
tions at the cut-off layer, at the density "bump" or 
after the fast mode reflection above the cs = va 
layer. In addition to these physical reasons, an- 
other nodes can be formed in the line-of-sight ve- 
locity projection due to the interference of the fast 
and slow magneto-acoustic modes. These type of 
nodes are more frequent at latitudes close to mag- 
netic equator where the amplitudes of the modes 
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become of the same order of magnitude. 

• In the atmosphere above logrs = the amphtude of 
the slow mode waves increases exponentiahy with 
height, in accordance with behavior of acoustic- 
gravity waves in a stratified atmosphere. Their 
phase behavior depends on the frequency and the 
field inclination. The amplitude and the phase of 
the fast mode waves vary only slowly with height 
in agreement with their nearly standing character. 

• In the strong-field case, when the whole simulation 
domain is in the field-dominated regime {Bq > 1 
kG), the low- frequency waves (below the cut-off) 
are evanescent in the atmosphere, except for the re- 
gions close to the magnetic equator, where the cut- 
off frequency is lowered due to their field-aligned 
propagation. The high-frequency waves (above the 
cut-off) are propagating at all magnetic latitudes 
in the strong-field case. The presence of the mode 
transformation changes this situation. In the weak- 
field case (^0 < 0.5 kG) both the high-frequency 
and low-frequency waves become partially trapped 
near the magnetic equator due to the refiections 
of the fast magneto-acoustic mode produced in the 
magnetically-dominated part of the atmosphere. 

• The temperature and density variations show an 
amplitude and phase jump around logrs = 0, where 
the properties of the atmosphere change in an 
abrupt way. The amplitudes of temperature and 
density variations are unusually high at these layers 
and can reach some 5-10% of the static values. The 
amplitudes of the magnetic field variations reach 
at most 3-4 G, decreasing rapidly with height and 
with the Bq field strength of the simulation. 

• The disc-integrated variations of the line-of-sight 
velocity and thermodynamic parameters mostly 
preserve signatures of the slow magneto-acoustic 
waves close to the magnetic poles for all angles be- 
tween the line-of-sight and and pulsation axis. The 
amplitudes of the variations increase exponentially 
with height. The node layer in velocity is main- 
tained around logrs = 0. The prominent increase 
in the amplitudes of variations of temperature and 
density around logrs = is also preserved in the 
disc-integrated data. 

Our conclusions regarding the existence of the fast and 
slow magneto-acoustic waves and their propagation prop- 
erties a re similar to t hose obtained in the theoretical 
work of IS^sa Cunhal (120081 ). The aim of our anal- 
ysis was to understand what mechanisms produce such 
diverse signatures of pulsations measured in the atmo- 
spheres of roAp stars. ISousa fc Cunhal (|2QQ8l ) addressed 
a different question, namely the energy losses of the high- 
frequency waves associated with the presence of magnetic 
field. Despite the different objectives, our modelling al- 
lows us to make some conclusions regarding the behav- 
ior of the high-frequency waves as well. In particular, 
we conclude that, if on their way to the surfaces the 
high-frequency fast magneto-acoustic waves (essentially 
acoustic in the interior where cs > va) travel through 
the region of equal Alfven and sound speeds where the 
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Fig. 15. — Example of the height-time variations of the vertical 
velocity for the global model with = 3 kG and Tq = 600 s at the 
magnetic latitude corresponding to ip = 62.7°. Yellow and white 
colors mean positive velocity, while dark red and blue colors mean 
negative velocity. The horizontal dotted lines correspond to the 
optical depths logrs = 1.5, 0, —2 and —4. Note the presence of a 
node around logrs = — 2 due to the fast and slow mode interference. 



mode transformation can occur, they can be partially 
transmitted as fast magnetic waves in the magnetically 
dominated atmosphere {cs < va)- The fast magnetic 
waves are then reflected back to the interior and thus, 
some part of their energy is lost at the surface. This way 
the high-frequency waves can become partially trapped. 
This mechanism acts more efficiently at high magnetic 
latitudes. The part of the energy of the acoustic-like fast 
modes in the interior transmitted as fast magnetic modes 
increases with increasing the inclination angle between 
the magnetic fleld and the direction of the wave prop- 
agation, thus becoming larger at the magnetic equator. 
At the poles where the fleld is close to vertical, the high- 
frequency fast acoustic-like modes from the interior are 
transformed completely into the slow acoustic-like modes 
in the atmosphere and no trapping occurs. 

T he difference b e tween our modelling and the one 
by ISousa Cunhal (|2008r ) is that we considered much 
smaller portion of the stellar interior. Thus, in our sim- 
ulations the wave transformation happens only for the 
weakest fleld strengths. Due to that, we have concluded 
that the high-frequency waves are trapped near the mag- 
netic equator only in the weakest flelds, below 0.5 kG. If 
we would have considered a larger portion of the interior 
of the star, the mode transformation could have affected 
waves for larger magnetic fleld strengths as well. 

The low-frequency waves with frequencies below the 
cut-off frequency are shown to be trapped at the mag- 
netic poles due to the cut-off effects. The low- frequency 
waves are also partially trapped near the magnetic equa- 
tor, since they are affected by the mode transformation 
in a similar way as high-frequency waves. However, it is 
possible that some part of their energy escapes at inter- 
mediate magnetic latitudes where the fleld inclination is 
already sufficient to lower down the cut-off frequency and 
the fast to fast mode transmissio n is not complete. Th ese 
conclusions are all in line with S ousa fc Cunhal (|2008l ). 

Our analysis allows us to link directly the theoreti- 
cally calculated wave perturbations in all quantities to 
the heights in the atmosphere, and, thus, to compare the 
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simulated and observed wave behavior. We stress that 
the conclusions coming from our analysis remain valid, 
despite we take only a small portion of the stellar interior. 
The relative distribution of the wave amplitudes at differ- 
ent frequencies may change due to the mode transforma- 
tion, but the general picture of the propagation, ampli- 
tude increase and the node formation will remain valid. 
For example, the consequences of the wave reflection 
(such as the standing wave formation or the downward 
wave propagation due to the mode transformation) can 
hardly be detected in observations since these processes 
happen below the visible surface in the our roAp star 
model. Only for the weakest field strengths {Bq < 0.5 
kG) these features can, in principle, be detected in the 
deepest observable atmospheric layers. 

Our disc-integrated data suggests that, in the assumed 
excitation model, the variations of all quantities in the 
upper atmosphere (except for the magnetic field) are 
dominated by slow ma gneto- acoustic niodes. This con- 
clusion is different from lSousa fc Cunhal (|2008f ) who con- 
sidered disc-integrated velocity from their simplified, an- 
alytical model for the case of the observer being pole- 
on in relation to the magnetic field axis. They suggest 
that the observed variations are the superposition of the 
acoustic running waves and magnetic standing waves and 
the phase depends on the relative contribution of both 
components. But for the magnetic and pulsation geome- 
try adopted in our study mostly the signatures of the slow 
acoustic-like waves remain in the disc-integrated data. 

One of the main goals of our work is to interpret the 
complex observational picture of the spectroscopic vari- 
ability of roAp stars. However, detailed analysis of ob- 
servations requires computation of the hydrodynamical 
models and spectrum synthesis for specific stars, taking 
into account their pulsation frequencies, magnetic field 
strengths, orientation of pulsation and rotation axes, and 
individual atmospheric structure. We plan to undertake 
such an in-depth investigation of selected roAp stars in 
the future studies. Here we limit ourselves to the qualita- 
tive comparison of our model predictions with the main 
observational features of roAp pulsations. 

The growth of the pulsational amplitude with height, 
observed in all roAp stars, is frequently considered in 
the context of a simple picture where the kinetic wave 
energy, pv'^ is approximately constant with height, so 
that pulsation velocity amplitude increases due to the ex- 
pone ntial decrease of d ensity in the stratified atmosphere 
(e.g. jKurtz et al.ll2007>) . Following these arguments, one 
can estimate an increase of velocity amplitude by a factor 
of 5-10 from Fe-peak to REE lines. This is not consis- 
tent with the observational studies, which often infer the 
REE to Fe-peak velocity amplitude ratios of 30-100. Our 
simulations demonstrate that the variations in the lay- 
ers probed by the lines of light and iron-peak elements 
are strongly influenced by the density inversion around 
logTs = 0. This structure, located at the base of the 
photosphere of all dwarf stars in the roAp temperature 
range, leads to the formation of the velocity node right in 
the region where lines with weak or undetectable pulsa- 
tional variation are formed. This node is also present in 
the disc-integrated data in the same layers, while higher 
in the atmosphere the amplitude of the line-of-sight ve- 
locity increases exponentially. Consequently, as the wave 
propagates outside the node region, the velocity ampli- 



tude increases by a factor of 30-50 in our simulations, in 
reasonable agreement with observations. 

The atmospheric structure around logrs = has 
equally important repercussions for the pulsational tem- 
perature variations. Our simulations predict promi- 
nent temperature changes in this region, which will 
very likely dominate the observed luminosity variation 
of roAp stars. Thus, our magneto- hydrodynamical mod- 
els offer a new perspective on the photometric variabil- 
ity of pulsating Ap stars. We suggest that its origin 
is not a smooth, weak temperature fluctuation, grad- 
ually changing with height almo st equivalently to Teff 
variation (Med upe fc Kurt jll998l ). but a high- amplitude, 
vertically-localized phenomenon. A complete under- 
standing of this temperature variation is crucial for re- 
lating spectroscopic and photometric time-resolved ob- 
servations of roAp stars (|Rvabchikova et al]l2007cl ) and 
for interpreting photom etric pulsations in different bands 
([Matthews et al.l [19961 ). The position and shape of the 
density inversion can be also an important factor de- 
termining the visibility of the roAp pulsations in the 
Hertzsprung- Russel l diagram. Excitation theories (e.g., 
iTheado et al.l l2009l ) predict roAp-type oscillations in 
stars much hotter than the empirical blue border of the 
roAp instability strip at Teff ^ 8100 K. It is possible that 
the absence of the density inversion at the right atmo- 
spheric height diminishes luminosity variations, making 
them undetectable for the photometric pulsation surveys 
of hotter Ap stars. 

The depth dependence of the phase of the velocity 
variations obtained in our simulations shows a broad 
agreement with observations. The pulsation phase in- 
creases with height, showing the pattern of the outward- 
running wave clearly seen in many roAp stars. The 
high-frequency models predict a substantial difference 
in the phase of REE lines formed at different layers in 
the upper atmosphere. The low- frequency models of- 
ten demonstrate a slow height variation of the pulsation 
phase, followed by a more rapid changes in the higher 
layers. This is reminiscent of the standing-like followed 
by the running-like wave behaviour reported in observa- 
tional studies. 

For a few roAp stars the radial velocity studies of dif- 
ferent REEs reveal a sudden drop of the pulsation am- 
plitude and a tt radian change in phase at high atmo- 
spheric layers. Occasionally, such high- lying node sur- 
faces in the vertical velocity appear in our simulations at 
high magnetic latitudes. These nodes are produced due 
to a superposition of the fast and slow magneto-acoustic 
modes that attain similar amplitudes sufficiently far from 
the pulsation poles. An example of the vertical veloc- 
ity variation at the magnetic latitude corresponding to 
iIj = 62.7° in the simulations with B^^ = 3 kG showing 
a high-lying node is given in Fig. [151 Though such kind 
of nodes are quite persistent at high latitudes, they are 
not retained in the integrated velocity data. It happens 
because the amplitude of pulsations is small far from the 
poles for the assumed mode geometry. However, we can 
speculated that the upper atmospheric nodes can become 
visible in the radial velocity data due to an inhomoge- 
neous distribution of chemical elements, located in spots 
close to the magnetic equator. Additionally, contribu- 
tion of these regions to the disk-integrated signal could 
be amplified by the magnetic distortion of the mode ge- 
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ometry (|SaiQll2QQ5f ). 

The non- detection of the pulsational variations of mag- 
netic field is in fine with the results of our numerical sim- 
ulations. The disk-integrated amplitudes of the fiuctua- 
tions of the field modulus and longitudinal field are well 
below the current detection limits attained in the time- 
resolved spectroscopic and spectropolarimetric observa- 
tions of roAp stars. Furthermore, our models show that 
the depth-dependence of the magnetic field variations is 
opposite to that of the radial velocity. Minute fiuctua- 
tions of magnetic field rapidly decrease above logrs = 0. 
Thus, an eventual detection of the pulsational varia- 
tion of magnetic field is far more likely using the lines 
formed in deep layers, characterized by very weak or un- 
detectable velocity oscillations. 

Finally, theoretical pulsation models computed in our 
study allow us to comment on the REE line profile varia- 
tion (LPV) of the REE lines forming in the upper atmo- 
spheres of roAp stars. This interesting topic was explored 
in several recent theoretical and observational studies. 
With the aim to expl ain the anomalous asymmetric LPV 
of strong REE lines, IShibahashi et al.l (2008) proposed a 
shock- wave model according to which the intrinsic pulsa- 
tion velocity amplitude at the REE line formation heights 
(logTs = —4 : —5) reaches values comparable to the 
sound speed. To retain an agreement with the observed 
substantially subsonic, sinusoidal radial velocity curves, 
Shibahashi et al. postulated that the REE line formation 
zone spans over a shock- wave train, so that the observed 
signal corresponds to the vertically-averaged pulsational 
variation. Our hydrodynamical models do not confirm 
this picture. Apart from a moderate non-linearity in the 
uppermost layers, the velocity variations are harmonic 



and their amplitudes are well below the sound speed. For 
a broad range of the magnetic inclinations {tp = — 80°) 
the REE line- forming zone above logrs = — 2 accommo- 
dates no more than one vertical pulsation wavelength. 
For very high field inclinations {tp > 80°) multiple nodes 
can appear in the vertical velocity but the contribution of 
the corresponding surface regions to the disk-integrated 
signal is negligible. 

An alternative phenome nological model of the p ulsa- 
tional LPV of REE lines ( Kochukhov et all 120071 pro- 
vides a closer match to the behaviour of our theoretical 
simulations. Kochukhov et al. showed that the observed 
LPV can be reproduced by a combination of the subsonic 
velocity changes and the line width variation. The lat- 
ter can be related to the temperature fiuctuation in the 
uppermost layers (Fig. [T3|). However, detailed spectrum 
synthesis is required to study in detail the LPV predicted 
by our models and to compare these predictions with the 
observed pulsational line behaviour. 

Detailed spectroscopic analysis and the comparison of 
the spectra of different chemical elements synthesized in 
our simulations with observations will be the aim of our 
future investigations. 
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